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We show that recursively generated Chebyshev expansions offer numerically efficient represen- 
tations for calculating zero-temperature spectral functions of one-dimensional lattice models using 
matrix product state (MPS) methods. The main features of this Chebychev matrix product state 
(CheMPS) approach are: (i) it achieves uniform resolution over the spectral function's entire spectral 
width; (ii) it can exploit the fact that the latter can be much smaller than the model's many-body 
bandwidth; (iii) it offers a well-controlled broadening scheme that allows finite-size effects to be 
either resolved or smeared out, as desired; (iv) it is based on using MPS tools to recursively calcu- 
late a succession of Chebychev vectors \t n ), (v) whose entanglement entropies were found to remain 
bounded with increasing recursion order n for all cases analyzed here; (vi) it distributes the total 
entanglement entropy that accumulates with increasing n over the set of Chebyshev vectors \t n ), 
which need not be combined into a single vector. In this way, the growth in entanglement entropy 
that usually limits density matrix renormalization group (DMRG) approaches is packaged into con- 
veniently manageable units. We present zero-temperature CheMPS results for the structure factor 
of spin-i antiferromagnetic Heisenberg chains and perform a detailed finite-size analysis. Making 
comparisons to three benchmark methods, we find that CheMPS (1) yields results comparable in 
quality to those of correction vector DMRG, at dramatically reduced numerical cost; (2) agrees well 
with Bethe Ansatz results for an infinite system, within the limitations expected for numerics on 
finite systems; (3) can also be applied in the time domain, where it has potential to serve as a viable 
alternative to time-dependent DMRG (in particular at finite temperatures). Finally, we present a 
detailed error analysis of CheMPS for the case of the noninteracting resonant level model. 

PACS numbers: 02.70.-c,75.10.Pq,75.40.Mg,78.20.Bh 



I. INTRODUCTION 

Consider a one-dimensional lattice model amenable to 
treatment by the density matrix renormalization group 
(DMRG),— ~— with Hamiltonian H, ground state |0) and 
ground state energy Eq. This paper is concerned with 
zero-temperature spectral functions of the form 

A bc (uj) = (0\BS(u}-H + E )C\0) , (1) 

which represents the Fourier transform J ^e lult G BC (t) of 
the correlator 

G BC (t) = <0|B(i)C(0)|0) . (2) 

One possible framework for calculating such spectral 
functions is to expand them in terms of Chebychev poly- 
nomials, as advocated in Ref. H. Such a Chebyshev ex- 
pansion offers precise and convenient control of the accu- 
racy and resolution with which a spectral function is to be 
computed. This is very useful, particularly when broad- 
ening the spectral function of a length- 1 system, which 
exhibits finite-size subpeaks with spacing wi ~ 1/1, in 
order to mimick that of an infinite system: if the lat- 
ter has structures (e.g. sharp or diverging peaks) which 
are not yet properly resolved at the scale ujl, the broad- 
ened version of the finite-size spectral function inevitably 
bears I-dependent errors in the vicinity of these struc- 
tures. Hence, when calculating the finite-size version of 
these structures for the length-L system, there is no need 



to achieve an accuracy beyond that of the expected In- 
dependent errors, and having convenient control of this 
accuracy can significantly reduce numerical costs. 

In this paper, we show that Chebyshev expansions 
offer numerically efficient representations for calculat- 
ing spectral functions using matrix product state (MPS) 
methods^riS, with numerical costs that compare fa- 
vorably to those of other established DMRG-based ap- 
proaches. In particular, the Chebychev MPS approach 
presented here, to be called CheMPS, allows the above- 
mentioned control of accuracy and resolution to be im- 
ported into the DMRG/MPS arena. 

The historically first approach for calculating spec- 
tral functions with DMRG is the continued-fraction 
expansion^ While this method requires only modest nu- 
merical resources, it is limited to low frequencies and it 
is difficult to produce reliable results with it in the case 
of continua (however, algorithmic improvements were re- 
ported recently-^). At present, the most accurate, but 
also most time-consuming approaches are: (i) the correc- 
tion vector (CV) method*^— and (ii) time-dependent 
DMRG (tDMRG) jLEilSriS in particular when combined 
with linear prediction techniques J^r— Since any new ap- 
proach must measure up to their standards, let us briefly 
summarize their key ideas, advantages and drawbacks. 

(i) To calculate A BC (cu) using the CV approach, it is 
expressed as 

A BC {u) = (0\B\C)„ , (3a) 
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in terms of the so-called correction vector 



\C) U — — lim — Im 



1 



H + E + iri 



e\o). (3b) 



The correction vector can be calculated (for finite broad- 
ening parameter rj) using cither conventional DMRG^2 ~— 
or variational matrix product state (MPS) methods^ A 
major advantage of this approach is that arbitrarily high 
spectral resolution can be achieved by reducing 77 and 
sampling enough frequency points. However, this comes 
at considerable numerical costs: first, a separate calcu- 
lation is required for every choice of oj (though in doing 
so, results for |C) w 's from previous frequencies can be in- 
corporated); and second, the calculation of \C) U involves 
an operator inversion problem that is numerically poorly 
conditioned, ever more so the smaller rj is. 

(ii) An alterative possibility is to use tDMRG to cal- 
culate the time-domain correlator G BC (t), Fourier trans- 
forming to the frequency domain only at the very end. 
To this end, one expresses 



G BL {t) =e^ at (Q\B\C) t 
in terms of the time-evolved state 
\C) t = e- iAt C\0) 



(4a) 



(4b) 



and uses tDMRG to calculate the latter. Two attractive 
features of this strategy are: first, it builds on an exten- 
sive body of algorithmic knowledge for efficiently calcu- 
lating time-evolution^ ' 16 ' 17 and second, a simple linear- 
prediction scheme^— can be used to extrapolate the 
time-dependence calculated for short and intermediate 
time scales to longer times, thereby improving the qual- 
ity of results at low frequency at hardly any additional 
numerical cost. However, obtaining reliable results over 
a sufficiently large time interval can, in itself, be nu- 
merically very expensive, since the time-evolution of the 
many-body state \C)t is accompanied by a strong growth 
in entanglement entropy. This unavoidably also implies 
a growth of tDMRG truncation errors. 

Note that in both of the schemes outlined above, signif- 
icant (often heroic) amounts of numerical resources are 
devoted to calculating a single state, \C) U for given u 
or \C)t for given t, as accurately as possible; the over- 
laps or expectation values of interest, namely (0|B|C) w 
for (0|j6|C)t, are only calculated at the end, in a single, 
final step, after \C) U or \C) t have been fully determined. 
Actually, these states are calculated so accurately that 
they would have been equally suitable for calculating any 
other quantity (correlator or matrix element) involving 
that state. In a sense, DMRG is asked to work harder 
than necessary: it is used to calculate a single state with 
"general-purpose accuracy" , whereas the accurate calcu- 
lation of a particular expectation value involving that 
state would have been sufficient. 

The main motivation for the present work is to at- 
tempt to reduce this calculational overhead by employing 
a representation of the spectral function that avoids the 
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Fig. 1. (a) Sketch of a spectral function whose spectral width 
Wa is much smaller than the many-body bandwidth W. Be- 
fore making a Chebychev expansion, we rescale the interval 
w G [0, W*\, with effective bandwidth W* = 2Wa, onto the 
interval oj' G [— W',W], shown in (b), with rescaled half- 
bandwidth W' = 1 — |et and a safety factor e t ~ 0.025. 



need for calculating a single state with such high accuracy 
and instead allows numerical resources to be focussed di- 
rectly on the calculation of the relevant expectation val- 
ues. This can be achieved by representing the spectral 
function via a Chebychev expansion ^ 2 ^ 25 whose coef- 
ficients, the so-called Chebyshev moments, can be cal- 
culated recursively using MPS tools. Below, we briefly 
summarize the structure and main features of such an ex- 
pansion, thereby providing both an introduction and an 
overview of the material developed in detail in the main 
part of this paper. 

The Chebychev polynomials T n (x) form an orthonor- 
mal set of polynomials on the interval x G [—1,1]. They 
are very well studied mathematically,— ~— and are widely 
used for function expansions since they have very favor- 
able convergence properties. As will be described in de- 
tail below, the spectral function can be represented ap- 
proximately by a so-called Chebychev expansion, which 
becomes exact for N — > 00, of the following form: 



Af{u) 



2W/W* 

W 1 - W' 2 



N-l 

50/iO + 2 X! 9nHnT n (uj') 
n=l 



(5) 



Here the Chebyshev moments fi„ = (0\B\t n ) are obtained 
from the Chebychev vectors \t„) = T„(H')C\0) , and the 
g n are known damping factors that influence broaden- 
ing effects. The primes indicate that the Hamiltonian 
H and frequency ui were expressed in terms of rescaled 
and shifted versions, H' and u>', in such a manner that 
an interval u> G [0, W*] that contains the entire spectral 
weight is mapped onto a rescaled band lu' G [— W , W] 
of halfwidth W < 1. 

This representation has several useful features: 

(i) It resolves the interval uj G [0, W*\ with a uniform res- 
olution of 0(W*/N). 

(ii) The range of frequencies over which the spectral func- 
tion has nonzero weight, say Wa (to be called its spectral 
width) is often significantly smaller than the many-body 
bandwidth of the Hamiltonian, say W, as depicted in 
Fig[U By choosing the effective bandwidth W* to be of 
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order Wa instead of W, huge gains in resolution are pos- 
sible. 

(iii) A well-controlled broadening scheme, encoded in the 
damping factors g n , is available that allows finite-size ef- 
fects to be either resolved or smeared out, as desired. 

(iv) The Chebyshev vectors \t n ) are calculated using a 
(numerically stable) recursion scheme, which exploits 
Chcbychcv recurrence relations to calculate |t„) from 
H'\t n -i) and \t n -2) (sec Eq. (j30f ). Thus, the expectation 
values from which the spectral function is constructed are 
built up in a series of recursive steps (see Eq. ([7]) below) 
instead of being calculated at the end in one final step. 

(v) The bond entropy of successive Chebyshev vectors 
\t n ) is found empirically to remain bounded with increas- 
ing recursion number n, thus the complexity of these vec- 
tors remains managablc up to arbitrarily large n. 

(vi) Finally, and from the perspective of numerical costs, 
most importantly: ChcMPS efficiently copes with the 
growth in bond entropy with increasing iteration num- 
ber that usually limits DMRG approaches. It does so by 
distributing this entropy over all \t n ), thereby packaging 
it into managable units (sec (v)). In particular, when 
constructing and using the states \t n ), one never needs 
to know more than three at a time (and after use may 
delete them from memory). Hence, it is not necessary to 
combine all information contained in all \t„) into a single 
MPS. 

Let us constrast this with the CV or tDMRG ap- 
proaches: imagine expanding the correction vector or 
time-evolved state in terms of the Chebyshev vectors \t n ) , 
i.e. expressing them as linear combinations of the form 



N-l 



IO^E c "l**>> \c) t 



N-l 



(0) 



respectively. (The coefficients C™ and C" are related 
by Fourier transformation.) Now, the CV or tDMRG 
approaches in effect attempt to accurately represent the 
entire linear combination using a single MPS. This en- 
devour is numerically very costly, since the entanglement 
entropy of this linear combination grows rapidly with N. 
The Chebychev approach avoids this problem by taking 
expectation values before performing the sum on n: 



N-l 



N-l 



(0|S|C\ ~ £ C>„ , (0\B\C) t ~ £ CtVn ■ (7) 

n=0 n=0 

Thus, the Chebychev expansion very conveniently orga- 
nizes the calculation into many separate, and hence nu- 
merically less costly, packages or subunits. 

Our paper is organized as follows. We introduce the 
Chebyshev expansion for spectral functions in Sec. Inland 
discuss its implementation using MPS including a new al- 
gorithm for performing a projection in energy in Sec. lIIII 
In Sec. IIVI we present CheMPS results for the structure 
factor of a spin-^ Heisenberg chains, perform a detailed 
analysis of finite-size effects (see Fig. [5]), and compare our 
results to CV, Bcthe Ansatz and tDMRG (see Figs. H M 



and respectively). In Sec. [V] we perform an exten- 
sive error analysis of the CheMPS approach using the 
quadratic resonant level model, and discuss some salient 
features of density matrix eigenspectra in Sec. I VII Sec- 
tion |VlT] summarizes our main conclusions, and Scc. lVIlJ 
presents a brief outlook towards possible future appli- 
cations, involving time dependence or finitc-tcmpcraturc 
correlators. An Appendix gives a detailed account of 
CheMPS results for the resonant level model used for 
the error analysis of Sec. [V] 



II. CHEBYSHEV EXPANSION OF A bc (oj) 

A. Chebyshev basics 

Let us start by briefly summarizing those properties of 
Chebyshev polynomials that will be needed below. We 
follow the notation of Ref. d, which gives an excellent 
general discussion of Chebyshev expansion techniques 
(though without mentioning possible DMRG/MPS ap- 
plications). 

Chebyshev polynomials of the first kind, T n (x), hence- 
forth simply called Chebyshev polynomials, are defined 
by the recurrence relations 



T n+ i(x) = 2xT n (x) - T n _i(x), 
T Q (x) = 1, Ti(x) = x . 



(8) 



They also satisfy the useful relation (for n > n') 

T n+n ,(x) = 2T n {x)T n ,{x) - T n _ n ,(x) . (9) 

Two useful explicit representations arc: 

T n (x) = cos [n arccos(cc)] = cosh [n arccosh(cc)] . (10) 

On the interval I = [—1,1] the Chebyshev polynomials 
constitute an orthogonal system of polynomials (over a 
weight function (irvl — x 2 )^ 1 ), in terms of which any 
piecewise smooth and continuous function f(x)\ x€l can 
be expanded. In fact, the T n {x) are optimally suited for 
this purpose, since they have the unique property (setting 
them apart from other systems of orthogonal polynomi- 
als) that on / their values are confined to |T„(x)| < 1, 
with all extremal values equal to 1 or — 1 . This is evident 
from the first equality in Eq. (p~0|) ; the second equality 
implies that for x (fc I, \T n (x)\ grows rapidly with in- 
creasing These properties are illustrated in Fig. [2] 

There are several ways of constructing Chebyshev ap- 
proximations for f(x)\ xeI (see Weisse et al.,— , Section 
II. A). The Chebyshev expansion that is practical for 
present purposes has the form 



/(*) 
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fj, n T n (x) 



(11) 
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(a) 




Fig. 2. (Color online) Chebyshev polynomials of the first kind, 
T n (x), for n up to 8. (a) All zeros and extrema of every T n (x) 
are located within the interval / = [—1,1], and all extremal 
values equal 1 or — 1. (b) Chebyshev polynomials |T n (a;)| for 
x e [—1.5, 1.5]. The |r„>o(x)| grow rapidly when \x\ increases 
beyond 1. 



where the Chebychev moments fj, n are given by 
l 

fi n = J dxf(x)T n (x) . (12) 
-l 

An approximate representation of order N is obtained 
for f{x) if only the first A" terms (i.e. n < N — 1) are re- 
tained. However, such a truncation in general introduces 
artificial oscillations, of period ~ 1/iV, called Gibbs os- 
cillations. These can be smoothened by employing cer- 
tain broadening kernels, which in effect rearrange the in- 
finite series (fTTj) before truncation. This leads to a re- 
constructed expansion of the form 



In(x) 



1 



ttVT 



jv-1 



.90^0 



+ 2 ^2 g n ^nT n (x) 



(13) 



which (for properly chosen kernels) converges uniformly: 
ma* ■ \f(x)-f N (x)\ N ^?0. (14) 

— Kx<l 



The reconstructed series ([13]) contains the same Cheby- 
shev moments \i n as Eq. (fT2"j) . but they are multiplied by 
damping factors g n , real numbers whose form is charac- 
teristic of the chosen kernel. Several choices have been 
proposed, which damp out Gibbs oscillations in some- 
what different ways (see Ref. H for details). We will 



mostly employ Jackson damping, given by 

(N - n + 1) cos ^ + sin ^ cot ^ 



N + 1 



(15) 



This is usually the best choice, since it guarantees an in- 
tegrated error of O(jj) for /jv(x). When used to approxi- 
mate a (5-function S(x—x) sitting at x € /, Jackson damp- 
ing yields a nearly Gaussian peak of width \/l — x 2 n/N. 
On one occasion we will also employ Lorentz damping, 



9n,X 



sinh [A 



1 



A 



)] 



sinh A 



(16) 



where A is a real parameter. Lorentz damping preserves 
analytical properties (causality) of Green's function and 
broadens a (5-function 6(x — x) into a peak whose shape, 
for the choice A = 4 used here (following Ref. , is nearly 
Lorentzian, of width \/l — x 2 X/N. 

To summarize: the order- N Chebychev-reconstruction 
/at(x) with Jackson or Lorentzian damping with A = 4 
yields a result that is very close to the broadened function 



dxK$ (x-x)f{x), 

'N .x 



(17) 



(X = J, L) with broadening kernels and widths given by 



-* 2 /(2(V 2 



V27T77' 

ffhf_ 

.2 + ^2 ' 



v'n.s = Vl-a; 2 ^ , (18s 



N.x — vl 



Vn 



N ' 



(18b) 



respectively. Thus, /jv(x) resolves the shape of f(x) with 
a resolution of 0(1/N). 

For purposes of illustration, Fig. [3]Ja) shows three 
Chebyshev reconstructions of a 5-f unction at x = 0: 
without damping, giving Gibbs oscillations; with Jack- 
son damping, yielding a near-Gaussian peak; and with 
Lorentz damping, yielding a near-Lorentzian peak. Fig- 
ure[3fb) shows a Jackson-damped Chebyshev reconstruc- 
tion of a comb of Gaussian peaks ^ a K?, {x — x a ), whose 
widths fj' a are all equal. It illustrates how increasing N 
reduces the amount of broadening until the original peak 
form is recovered for sufficiently large N. It also shows 
that the broadened peak widths depend on the peak po- 
sitions, reflecting the fact that convolving a Gaussian of 
width fj' a with a near-Gaussian of width rj' N Sq (Eq. (fT7|) ) 
produces a near-Gaussian of width 



v'l. 



(19) 



To evaluate the T n {x) that occur in Eq. (P~5|) . we use the 
first equality of Eq. (fTO]) . Although numerically more ef- 
ficient methods exist for this purpose,— their use becomes 
advisable only for expansion orders much larger than the 
A" < O(10 3 ) that we will need in this work. 
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Fig. 3. Three Chebyshev reconstructions of S(x), with N = 
50: the undamped case (g n = 1) yields Gibbs oscillations 
(central peak has height 5so(0) = 16.23); Jackson damping 
(S J ) mimicks a Gaussian peak K J (x) of width tt/N; Lorentz 
damping (5 ) for A = 4 mimicks a Lorentzian peak K L {x) of 
width X/N. Inset: Jackson and Lorentz damping factors, 
and <7,^a=4! respectively, plotted for N = 50. Both decrease 
monotonically from 1 to 0, but in somewhat different ways, 
(b) Jackson-damped reconstruction of a comb of normalized 
Gaussians (dashed line), all of width fj' = 0.02, for three val- 
ues of N (solid lines) . The ^-dependence of the peak heights is 
given by [2Ti(fj' 2 + r]'^ x )]~ 1/2 (dash-dotted line), see Eq. (fl9|> . 



slightly smaller than the requisite [—1,1] by choosing the 
rescaled half-bandwidth W to be smaller than 1, with a 
safety factor— of e t ~ 0.025. To be explicit, wc define 



& , = H_Eo _ w , 



(21a) 
(21b) 



where H' has ground state energy E' = —W. Then we 
express the spectral function ([1]) as 



A bc {lj) = -(0\B6{u' - H')C\0) 
a 



(22) 



(with u)' = lo'{lo) and H' given by Eqs. (f2TjO . which by 
construction has no weight for lj' ^ [— W', W']. 

One possible choice for W* is to equate it to the width 
of the many-body spectrum of H, given by W = E miac — 
Eq. When using DMRG, Eq is usually already known 
from calculating the ground state |0) of H, and E max can 
be found, e. g. , by calculating^ the ground state of — H 
(reduced DMRG accuracy relative to usual ground state 
calculations is sufficient, since only E max is of interest 
here.) 

A disadvantage of the choice W* — W is that the many- 
body bandwidth W typically is large (it scales with sys- 
tem size), whereas optimal spectral resolution requires 
W* to be as small as possible: since an N-th order Cheby- 
shev expansion yields a resolution of 0(1/N) on the inter- 
val [—1, 1], its resolution on the original interval [0, W*\ 
will be 0(W*/N), which evidently becomes better the 
smaller W*. If B and C are single-particle operators, the 
spectral width Wa of A bc (lj) is independent of system 
size and hence much smaller than the many-body band- 
width W . In this case, it is advisable to choose W, to be 
of similar order (though still larger) than Wa. We will 
choose W* = 2Wa, which is typically <C W, as illustrated 
in Fig. [TJ 



B. Rescaling of uj and H 



C. Chebychev expansion in frequency domain 



To construct a Chebyshev expansion of the spectral 
function A bc (uj) of Eq. ©, we need to rescale and shifts 
the Hamiltonian H i— > H' and the frequency lj t— > uj' in 
such a way that the spectral range of A(uj), i.e. the inter- 
val [0, Wj\\ within which it has nonzero weight, is mapped 
into the interval [—1, 1]. Rescaled, dimensionlcss energies 
and frequencies will always carry primes. As safeguards 
against "leakage" beyond [—1,1] due to numerical inac- 
curacies, we choose the linear map (see Fig. [I} 



e [0, W*\ h-> u' g [-W, W] , W = l 



(20) 



which entails two precautionary measures: first, the u- 
interval is taken to be larger than the requisite [0, Wa] by 
choosing the effective bandwidth W* to be larger than the 
spectral width Wa; second, the w'-interval is taken to be 



To expand the (5-function in Eq. (|22|) in Chebyshev 
polynomials, we use f{x) = 8{x — H') with x = lj' 
in Eq. (fT2|) . and obtain from Eq. (fT3]) a reconstructed 
Chebyshev operator expansion of the form: 



6 n (lj'-H') = 



1 



ttVI 



9o 



N-l 

! E 

n=l 



g n T n (H')T n (Lj') 



(23) 

Inserting this into Eq. (|22|) for A bc (oj) yields the Cheby- 
shev expansion ([5]), with Chebyshev moments given by 



(0\BT n (H')C\0) 



(24) 



Thus fi n is a ground state expectation value of an n-th or- 
der polynomial in H\ whose construction might a priori 
appear to become increasingly daunting as n increases. 
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Fortunately, this challenge can be dealt with recursively, 
by expressing the moments as 

= (0\B\t n ) , \t n ) = T n (H')C\0), (25) 

and calculating the Chebyshev vectors \t n ) by exploiting 
the Chebyshev recurrence relations ([5]). The details of 
this recursive scheme will be discussed in Section Hill 



D. Chebychev expansion in time domain 

The Chebyshev expansion can also be employed for 
studying time evolution in general, and the correlator 
G sc (t) in particular. To this end, we express the time- 
evolution operator as 



U(t) 



-iHt 



fa' e -i[*(.»'+w')+£io]t 5 ( u > _ h') , (26) 



and insert Eq. (|23[) (without damping, g n = 1) into the 
latter. This yields^! 



u N {t) 



Cn{t) 



,-i(E +aW')t 



N-l 



c (t) + 2^2T n (H')c n (t) 

- w r„(a/) 



n=l 



irVT 
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-duj 1 = (~i) n J n (at) 



(27a) 
(27b) 



Here J n {at) is the Bessel function of the first kind of order 
n. It decays very rapidly with n once n > at. Hence, 
an expansion of given order N gives an essentially exact 
representation of U(t) for times up to i max < — , while 
c N-i(t) provides an estimate of the error. 

Inserting Eqs. (gTJ) into Eqs. Q for G Be (t) we find 



Gfit) 



-iaW t 



N-l 



/'0 



J {at) + 2 ^H)>« J n {at) 



(28) 

where the Chebyshev moments [i n are again given by 
Eq. (|25jl . Thus the Chebyshev expansions of G BC (t) and 
A (cj) are governed by the same set of moments (j, n , as 
is to be expected for functions linked by Fourier trans- 
formation. 



III. MPS EVALUATION OF THE CHEBYSHEV 
MOMENTS fi n 

We now present a recursive scheme for calculat- 
ing the Chebyshev moments /i„. The manipulations 
described below were implemented using MPS-based 
methods,—^— which are very convenient for constructing 
the states of interest, while matrix-product operators^ 
(MPOs) simplify the implementation of the shift- and 
rescaling transformation Eq. (|21b() of the Hamiltonian. 



A. Recurrence fitting 

To initialize the Chebyshev expansion, we calculate 
ground state |0) and ground state energy E of H, make 
a specific choice for W* and W , and construct H ' accord- 
ing to Eq. (|21bp . Then comes the main task, namely the 
recursive calculation of the moments /i n . This is done 
starting from 

\t Q )=C\0), \t 1 )=H'\t ), (29) 

and using the recurrence relation (obtained from Eq. ([8])) 

\t n ) = 2H'\t n ^) - |t„_ 2 ) . (30) 

Eq. ([30)) can be implemented using the so-called compres- 
sion or fitting procedure^ (see^, Sec. 4.5.2 for details). It 
finds an MPS representation for \t n ), at minimal loss of 
information for given MPS dimension to, by variationally 
minimizing the fitting error 



A fit = \t n )- (2F'|t n _i)-|t n _s 



(31) 



We will call this procedure recurrence fitting. In practice, 
the variational minimization proceeds via a sequence of 
fitting sweeps back and forth along the chain. These 
are continued until the state being optimized becomes 
stationary, in the sense that the overlap 



1 



(t n | t n 



nun moil 



(32) 



between the states \t n ) and \t' n ) before and after one fit- 
ting sweep, drops below a specified fitting convergence 
threshold (typically in the range 10~ 6 to 10~ 8 ). The 
maximum expansion order for which \t n ) is obtained us- 
ing recurrence fitting will be denoted by -/V max . 

The MPS dimension to needed to achieve accurate re- 
currence fitting turns out to be surprisingly small (see 
Scc.[V]for a detailed analysis). For example, to = 32 suf- 
ficed for the antifcrromagnctic Heiscnbcrg chain of length 
L = 100 discussed in Section [TV] The reason for this re- 
markable and eminently useful feature lies in the fact that 
the Chebychev recurrence relations (|30| contain only two 
terms on the right-hand side, whose addition requires 
only modest computational effort. In contrast, CV or 
tDMRG typically require much larger to, since they at- 
tempt to represent the sum of many states, see Eq. (|6|). 
in terms of a single MPS. 

For the special but common case that B — C\ Eq. © 
yields a relation between different moments, 



(33) 



This can be used to effectively double the order of the 
expansion to 2A^ ma x without calculating any additional 
Chebyshev vectors, by setting n 1 = n — 1 or n: 



M2n-1 — 2(t n — /il, 

M2ri = 2(t„ \t n ) — Hq. 



(34) 
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We use tildes to distinguish /i„-momcnts calculated in 
this manner from the /i n -moments obtained via Eq. (|25jl . 
Although they should nominally be identical, in numer- 
ical practice /Et n -moments are less accurate (by up to a 
factor of 5 in Fig. [9jc) below), since they depend on 
two Chcbyshcv vectors, whereas /i„-moments depend on 
only one. Our Chebyshev reconstructions thus gener- 
ally employ the /i rl -moments, and unless stated other- 
wise, /i„-moments are used only for results requiring 



B. Energy truncation 

We have argued above that in order to optimize spec- 
tral resolution, it may be desirable to choose the effec- 
tive bandwidth W, to be smaller than the full many-body 
bandwidth W. If this is done, however, it is essential to 
include an additional energy truncation step into the re- 
cursion procedure, to ensure that each \t n ) remains free 
from "high-energy" components, i.e. ff'-eigenstates with 
cigcnencrgics E' k > 1, which fall outside the range [—1,1] 
that is admissable for arguments of Chebyshev polynomi- 
als. If < W, numerical noise causes the state \t n ) to 
contain such high-energy contributions in spite of the pre- 
cautionary measures described after Eq. because the 
application of H' to |t„_i) in Eq. ([30| entails a DMRG 
truncation step, which is not performed in the eigenba- 
sis of H' . If such high-energy components were fed into 
subsequent recursion steps, the norms of successive 

Chebyshev vectors would diverge rapidly (as would the 
resulting moments /i„), because this effectively amounts 
to evaluating Chebyshev polynomials T n (x) for \x\ > 1, 
where |T„| > 1 (sec Fig.^b)). 

As a consequence, after obtaining a new state \t n ) from 
Eq. (|30p . we take the precautionary measure of projecting 
out any high-energy components that it might contain, 
before proceeding to the next |£ n +i). This can be done 
by performing several energy truncation sweeps. During 
an energy truncation sweep, we focus on one site at a 
time, perform an energy truncation in a local Krylov ba- 
sis constructed for that site, and then move on to the next 
site. Shifting the current site is accomplished by stan- 
dard MPS means, without any truncation, as a DMRG 
truncation would counteract the energy truncation. (As 
a consequence, an energy truncation in terms of two-site 
sweeps has not been implemented.) 

The truncation must take place in the energy eigen- 
basis of the Hamiltonian H' . Of course, its complete 
cigenbasis is not accessible, thus we build a Krylov sub- 
space of dimension cLk within the effective Hilbert space 
at every site. Alternatively, energy truncation can also be 
performed in the bond representation \ip) = E>i r \lk)\rk) . 
In this Krylov subspace, the effective Hamiltonian H' K 
of dimension dx can be fully diagonalized and so we can 
construct a projection operator to project out all eigen- 
states with energy bigger than some energy truncation 
threshold ep . The choice of this threshold depends on the 



choice of M4- We have found the combination W* = 2Wa 
and ep = 1.0 to work well (but other choices, involving, 
e.g. smaller W* and larger ep would be possible, too.) 

In the following, we describe the procedure just out- 
lined in more detail for a single site, using standard MPS 
nomenclature. Let the effective local Hilbert space for 
this site be spanned by the left, local and right basis vec- 
tors \l), | cr) and |r), and expand the Chebyshev vector 
\tp) = \t n ) in this basis: 



(35) 



lar 



To construct a projection operator P that projects out 
the high energy components for this site, \ip) P\ip), 
one may proceed as follows: 

First, build a Krylov subspace of dimension dx within 
span{|Z) \o) \r}} and calculate the matrix elements of H' 
within it (no truncation necessary): 

\i) = {H'f- 1 \iP), i = l...d K , (36a) 

\i) i— > \i) orthonormalize via Gram-Schmidt, 

(36b) 

{H' K ) tj = (i\H' K \j), H>K eC d K *d K (36c) 

Next, fully diagonalizc H' K to obtain all eigenenergies e' a 
and eigenvectors \e a ) 



&H K U = Y^\e a )e' a (e a \. (37) 

a=l 

Then construct the projection operator 

p = i- e i e «>< e «i ( 38 ) 

a:£^>£_P 

for a certain energy threshold ep and apply it: 

\1>)»P\il>). (39) 

Performing this procedure once for every site of the chain 
constitutes a truncation sweep. The state obtained after 
several truncation sweeps, say \t n )tr, is stripped from the 
unwanted high-energy components of |t n ), as well as pos- 
sible within a Krylov approximation. After fitting and 
truncation have been completed, the resulting (unnor- 
malized) state |t n )tr is renamed \t n ), used for calculating 
/!„, and fed into the next recursion step. 

To quantify the effects of energy truncation, we con- 
sider two measures of how much \t n ) changes during trun- 
cation. First, for a given truncation sweep, we define the 
average truncated weight per site (averaged over all sites) 
by 



A rSWCCp 



(40) 



k a:e' >ep 
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parameter 


recommended value 


description 


task 


w* 


2Wa (or W) 


effective bandwidth with (or without) energy truncation 


1 Col^CLilll ti Ul J. J. 


et 


0.025 


safety offset in rescaled half-bandwidth: W' = 1 — ^et 


m 
A c 


10" 6 . . . HT 8 


MPS dimension 
fitting convergence threshold 


1 cL Hi 1 cllL-c 11 L Llllji 


(Ik 


30 


Krylov subspace dimension 




n s 


10 


number of sweeps 


energy truncation 


ep 


1.0 


energy truncation threshold (in rescaled units) 




N 

9n 


depends on system size 

gi 


order of expansion, broadening 
choice of damping factors 


spectral reconstruction 



TABLE I. List of CheMPS parameters that control various algorithmic tasks, influencing their numerical costs and the quality 
of results. Since the tasks "recurrence fitting" and "energy truncation" are carried out at every recursion step, the importance 
of the corresponding parameters is self-evident. However, W* and e t , which determine the rescaled Hamiltonian H' , turn out 
to have a high impact on the results, too, as the quality of the energy truncation sweeps strongly depends on H'. For ep = 1, 
the choice of taking W* to be twice the spectral bandwidth Wa (or equal to the many-body bandwidth W) was found to work 
well with (or without) energy truncation, respectively. N and g„ do not affect the calculated moments of the expansion, but 
control the broadening of the reconstructed spectral function. 



where |e^) are the vectors constituting the projector of 
Eq. (|38]l at site k. Second, we define the truncation- 
induced state change by 



At, = |||Utr- IUII 2 • (41) 



It measures changes in the state due to the intended trun- 
cation of high energy weight, but also due to unavoidable 
numerical errors. In our experience, neither of the trun- 
cation measures N™ eep and A tr show clear signs of decay 
when increasing the number of truncation sweeps, say ns 
(see Fig. HOT c) below). This reflects the fact that energy 
trunctation has the status of a precautionary measure, 
not a variational procedure, and implies that there is no 
dynamic criterion when to stop truncation sweeping. As 
a consequence, one has to analyze how the accuracy of 
the results depends on ns and optimize the latter accord- 
ingly. This will be described in Sec. IV Bl below. 

The numerical costs for energy truncation are as 
follows: The cost for the steps in Eqs. (|36|) are 
O (d 2 K m? d 2 D h) i where m is the MPS dimension, d the 
size of the local site basis and Dh the matrix product 
operator dimension of H'. The diagonalization of H' K 
is of 0(d K ) where dx is theoretically bounded by m 2 d. 
In our experience, the purpose of the energy truncation, 
which is solely to eliminate high-energy contributions, is 
well accomplished already for a relatively small Krylov 
subspace dimension of dx = 30 <C m 2 d. 

An overview of all the parameters relevant for CheMPS 
is given in Table [J Where applicable, it also lists the 
values that we found to be optimal. A detailed error 
analysis, tracing the effects of various choices for these 
parameters, will be presented in Sec. fVl 



IV. RESULTS: HEISENBERG 
ANTIFERROMAGNET 

To illustrate the capabilities and power of the pro- 
posed CheMPS approach, this section presents results 
for the spin structure factor of a one-dimensional spin- 
i Hciscnberg antiferromagnct (HAFM) and compares 
them against results obtained from CV and tDMRG ap- 
proaches. 



A. Spin structure factor 

We study the spin-i HAFM for a lattice of length L 

L-l 

-Hhafm = J 2^ ^ ' ' (42) 
i=i 

where Sj denotes the spin operator at site j. We 
choose J = 1 as unit of energy throughout this section. 
This model exhibits SU(2) symmetry, which has been 
exploited^ in our calculations , accordingly all MPS di- 
mensions noted for the HAFM are to be understood as 
number of SU(2) (representative) states being kept. To 
account for the open boundary conditions we define spin 
wave operators as: 

*-V^i>£i*- (43) 

The spin structure factor (spectral function) we are in- 
terested in is given by 

S(k,u) =A S *- Sk (u)- (44) 
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It is known from exact solutions^— that the dominant 
part of the spin structure factor stems from two-spinon 
contributions, bounded from below and above by 



u>i = — |sinfc| and ui2 = n 



Moreover, for an infinite system S(k, ui) is knowi 
diverge as 



(45) 



S(k,u) 
S(tt,uj) 



LOl\ 



uj- 1 v/ln(l Jw) 



■i^ln[l/(u-uji)] ,for/c^7r, (46a) 

(46b) 



as oj approaches the lower threshold ui\ from above. This 
divergence reflects the tendency towards staggered spin 
order of the ground state of the Hciscnberg antifcrro- 
magnct. It poses a severe challenge for numerics, which 
always deals with systems of finite size, and hence will 
never yield a true divergence. Instead, the divergence 
will be cut off at u> — u)\ ~ 1/L, yielding a peak of finite 
height, 



maxS'(A;,a;) 
max S(tt, uj) 



[L]nL] 
L[\nLY 



for k 7^ 7r, 



(47a) 
(47b) 



Thus, the best one can hope to achieve with numerics is 
to capture the nature of the divergence as u approaches 
Ui before it is cut off by finite size, or the scaling of the 
peak height with system size. 

Eq. (|45p gives a good guide for choosing W* . We found 
the choices W* = 6.3 ~ 2ir and et = 0.025 to work well 
for all k and have used them for all figures ([4] to [6|) of 
this section. As consistency checks, we verified that the 
resulting S(k,uj) is essentially independent of W*, and 
that it agrees with a calculation that included the full 
many-body bandwidth (W* = W). 

To have an accurate starting point for all calculations, 
we throughout used a ground state obtained by standard 
DMRG with MPS dimension m = 512. From expansion 
order n = 1 onwards, it turned out to be sufficient to 
represent all Chcbyshcv vectors \t n ) using a surprisingly 
small MPS dimension of m = 32, or m = 64 for some 
results involving very large iteration number, as indicated 
in every figure. (In retrospect, this implies that for the 
ground state, too, a much smaller m would have sufficed.) 
We have verified that the structure factor S(k,uj) is well 
converged w.r.t. m nevertheless. Detailed evidence for 
this claim will be presented below. However, already 
at this stage it is worth remarking that the ability of 
CheMPS to get good results with comparatively small m- 
values is perhaps the single most striking conclusion of 
our work. This will be discussed in detail below. 



B. Comparison to CV 

We begin our discussion of CheMPS results by compar- 
ing them to those of CV calculations, which are known 



3 
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Ch L :N=111 

Ch L : N=222 

CV:r,=0.1 + 
CV: n=0.05 X 

L=100, 1=4.0, a=3.19 
m cv =1000, m Ch =32 
W*=6.3, e t =0.025 



co/ji 



Fig. 4. (Color online) Comparison of CheMPS vs. correc- 
tion vector (CV) calculations of S(k = n/2,oj) for a Heisen- 
berg chain: Lines show Chebyshev results reconstructed for 
N = 118 and 236 using Lorentz damping, Eq. (|16p . with 
A = 4.0; symbols show CV results, obtained using broad- 
ening parameters of n — 0.1 and 0.05. We expect and in- 
deed find good agreement between lines and symbols, since 
both Lorentz damping and the correction vector method in 
effect broaden the spectral function by Lorentzians, whose 
widths we equated by choosing aX/N = n (with a = 3.19), 
see Eq. (|4"8")l , and also Fig. [3] Since mch "C "lev, the nu- 
merical cost of obtaining an entire curve via Chebychev is 
dramatically cheaper than calculating a single point via CV, 
as discussed in the text. 



to be very accurate, though also computationally expen- 
sive. The CV method involves a broadening parame- 
ter rj and broadens (^-functions into Lorentzian peaks of 
width rj. This can be mimicked with CheMPS by us- 
ing Lorentz damping (with A = 4.0), since this also pro- 
duces Lorentzian broadening, representing a (5-function 
8{ui' — uj') by a near-Lorentzian peak, albeit with a 
frequency-dependent width, r)' N Q , = VT — uj 12 X/N (see 
Sec. ILTAl and Fig.©. To compare CheMPS results with 
CV results at given 77, we thus identify T] = ar]' N -,, where 
a is the scaling factor from Eqs. (|2"Tj) and u>' is taken to 
be the rescaled and shifted version of the frequency w max 
at which the peak reaches its maximum. Thus we set the 
expansion order used for reconstruction to 



4a 



N=— v 7 ! - KWa - VF') 2 



(48) 



Figure 3] shows such a comparison for the structure fac- 
tor S(tt/2,uj) of a L = 100 Hciscnberg chain. Wc 
used two choices of ?y that are large enough to avoid 
finite-size effects, namely 77 = 0.1 and 0.05, and set 
Wmax = tt/2 (cf. wi of Eq. (|4T)])). We used MPS dimen- 
sions of m cv = 1000 or m ch = 32 for CV or CheMPS 
calculations, respectively. (Our choice for mcv aimed for 
achieving highly accurate CV results; for 77 = 0.05 this 
required mcv = 1000, but for 7/ = 0.1, a slightly smaller 
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value for tocv would have sufficed.) We find excellent 
agreement between the two approaches without adjust- 
ing any free parameter, since TV is fixed by Eq. (jJSJ). For 
example, for r\ = 0.05, N = 255, the relative error is less 
than 3 % for all ui. 

Since this level of agreement is obtained using mch "C 
mcv, we conclude that ChcMPS with Lorentz damping 
gives results whose accuracy is comparable to those of 
CV, at dramatically reduced numerical cost. Indeed, for 
i] = 0.05 the calculation of the entire CheMPS spectral 
function was 25 times faster than that of a single CV 
data point. 



C. Finite-size effects 

Let us now analyse the role of finite system size. To this 
end, it is of course important to understand broadening 
effects in detail. The fact that ChcMPS offers simple 
and systematic control of broadening via the choice of 
the expansion order (and damping factors), as will be 
illustrated below, is very convenient and may regarded 
as one of its main advantages. 

Figure [5ja) shows CheMPS results for the spin struc- 
ture factors S(k,u)) of four different momenta fc, calcu- 
lated for L = 100 using Jackson damping. They were 
reconstructed using the largest expansion order, say Nl, 
that does not yet resolve finite-size effects, a choice that 
will be called optimal broadening. Each curve shows a 
dominant peak, and we are interested in finding its in- 
trinsic shape S°°(k,uj) in the continuum limit of an in- 
finitely long chain (L — > oo). Thus, the following general 
question arises: under what conditions will a spectrum 
calculated for finite system size L and reconstructed with 
finite expansion order TV, say S^j(u>), correctly reproduce 
the desired continuum spectrum S , °°(w)? The general an- 
swer, of course, is that the optimally broadened spectrum 
should have converged as function of L, i.e. the shape 
of Sj^ L (uj) should not change upon increasing L. How- 
ever, for a spectrum with an intrinsic divergence, such as 
Eq. (|4"6")h the peak's height will never saturate with L; at 
best one can hope to observe L-convergence of the shape 
of its tail, and the proper scaling of its height (Eq. (07])). 

To illustrate the nature of finite-size effects and the 
role of TV in revealing or hiding them, Fig. [5jb) shows 
S(ir, ui) for L = 100 and several values of TV, both smaller 
and larger than TVl. As N is increased and the effective 
broadening ij N ~ (D(W*/N) decreases, the main peak 
of the initially very broad and smooth spectral function 
becomes sharper. Optimal broadening in Fig. [HJb) cor- 
responds to Nl — 70, beyond which additional "wiggles" 
emerge. These develop, with beautifully uniform resolu- 
tion, into dominant subpeaks as TV is increased further. 
The discrete subpeaks reflect the quantized energies of 
spin-wave excitations in a finite system. With sufficiently 
high resolution (TV = 999 in Fig. [UJb)) numerous addi- 
tional minor subpeaks emerge, but their weight is very 
small compared to that of the dominant subpeaks. This 



fact is important, since it implies that the structure fac- 
tor of a finite-size system is exhausted almost fully by 
the set of dominant subpeaks, with very small intrinsic 
widths. 

We have checked that there are 0(L) dominant sub- 
peaks within the spectral bandwidth of S{k,u). Corre- 
spondingly, the average spacing between dominant sub- 
peaks, to be called the finite-size energy scale lol, is pro- 
portional to \ (Fig. [Hcd)). The weight of each subpeak 
decreases similarly, ensuring that the total weight in a 
given frequency interval converges as L — > oo. The in- 
verse subpeak spacing H/ujl corresponds to the Heisen- 
berg time, i.e. the time within which a spin wave packet 
propagates the length of the system. 

Figures E)(e,f) illustrate two slightly different broaden- 
ing strategies. In Fig. [He), L is increased for fixed TV: 
the distinct subpeaks increasingly overlap, resulting in a 
smooth spectral function once wj, drops below tjn. In 
Fig-EJf) optimal broadening is used (t]n just larger than 
lul- now no subpeaks are visible, and the L-cvolution of 
the main peak is revealed with better resolution. 

In both Figs. [5](e) and[5{f), the peak height shows no 
indications of converging with increasing L. (The same 
is true for the data of Fig. [SJa).) This reflects the intrin- 
sic divergence of the peak height expected from Eq. (|46jl . 
Figures Efa) and &b) contain a quantitative analysis of 
this divergence, for S(tv,uj) and S(ir/2,uj), respectively. 
The shape of the divergencies for an infinite system are 
shown by the thick solid lines, representing exact Bethe 
Ansatz results from Ref. HH. Thin dashed lines show re- 
sults from tDMRG from Ref. H for L = 100, and thin 
solid lines CheMPS results for several system sizes be- 
tween L = 50 and 300. For CheMPS spectral recon- 
struction, we determined the expansion order TV300 that 
ensures optimal broadening for L = 300, and used a 
fixed ratio of N/L = TV300/3OO for all curves (namely 
0.42 or 0.67 for Figs. E[a,b) , respectively). CheMPS (for 
L = 300) and tDMRG reproduce the peak's tail and 
flank well, but clearly and expectedly are unable to pro- 
duce a true divergence at the lower threshold frequency. 
Nevertheless, the insets show that the manner in which 
the CheMPS peak heights increase with L is indeed con- 
sistent with Eq. (|46|) . (For the limited range of avail- 
able system sizes, however, a reliable distinction between 
Lfh^L)] 1 / 2 , [Lln(L)] 1 / 2 or L behavior is not possible.) 

It is also possible to determine the lower threshold 
frequency lo\ rather accurately from the CheMPS re- 
sults by doing an 1/L extrapolation. We illustrate this 
in Fig. EDJb) by extrapolating the frequencies at which 
S(it/2,l>j) = 0.1 (triangles). Since the data exhibit a 
slight curvature when plotted against 1/L (see lower in- 
set of Fig. Efb)), they were fitted using a second order 
polynomial in 1/L. Extrapolating the fit to 1/L = 
yields ui = 0.4967T (marked by a square), in good agree- 
ment with the prediction wi = n/2 from Eq. (1451) . 
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Fig. 5. (Color online) Spin structure factors for a Heisenberg chain, reconstructed using Jackson damping, (a) S(k, oj) for 
four choices of momentum k, for a chain of length L = 100. Each spectrum was reconstructed using optimal broadening, i.e. 
by choosing the largest expansion order, say Nl (indicated by superscripts), that does not yet resolve finite-size effects, (b-f) 
Finite-size analysis of S(tt,ui). (b) To determine Nl for given L (here 100), several different expansion orders N are considered. 
Increasing N reduces the effective broadening r]N — 0(W*/N) until finite-size subpeaks appear for N > Nl (here Nl ~ 70, bold 
red curve), (c) Evolution of the finite-size structure with L, revealed by fixing iV large enough (here = 499) to resolve the first 
few dominant subpeaks of all curves. There are L dominant subpeaks (not all shown here) within the spectral bandwidth, with 
average spacing ujl ~ 1/L. (d) Same as (c), but plotted on a semi-log scale, and with somewhat smaller N (here = 180), chosen 
to be somewhat larger than the optimal broadening Nl for the largest L (here N300 = 125). As L increases and uil decreases, 
the subpeaks coalesce toward the intrinsic lineshape S°°(k,Ld). (e) When L is increased at fixed N (here 70), finite-size effects 
disappear once ujl drops below the effective broadening r]N, resulting in a smooth spectral function, (f) In contrast, when L 
is increased while using optimal broadening, N = Nl (i.e. wn just above ujl), none of the curves show finite-size effects, and 
the resulting main peak is sharper than in (e). In both (e) and (f), the peak height shows no indications of converging with L, 
reflecting the fact that the true peak shape involves an w _1 [lri w] -1 ' 2 divergence. Moreover, the CheMPS curves in (f) show 
signs of overbroadening when compared to the exact Bethe Ansatz result (dashed), from Ref . |38|. 



D. Discrete representation of spectral function 

In both Fig. |5jf) and Fig. [6j the right flank of the 
peak still bears signatures of overbroadening: the curve 
for a given L lies above those for larger L (before bend- 
ing over towards its peak), and all curves lie significantly 
above the exact Bethe Ansatz curve (dashed line). One 
way of reducing this broadening would be to simply in- 
crease L, but this is numerically costly. Clearly, alter- 
native strategies for reducing finite-size effects would be 
desirable. One such scheme, involving linear prediction 
in the time domain, will be discussed in the next subsec- 
tion. Here we present another, which exploits the ability 



of CheMPS to accurately resolve finite-size peaks. 

The origin of overbroadening is clear: when neighbor- 
ing subpeaks are broadened enough to overlap, weight is 
inevitably transfcred from large peaks to smaller peaks. 
This effect is negligible only in the limit L — > 00, where 
the subpeak spacing becomes negligible. To avoid over- 
broadening for a finite- L system, one thus has to analyse 
spectra for which N is large enough that subpeaks do not 
overlap significantly, such as that shown in Fig.[5jb). 

To be concrete, let us represent the true, discrete spec- 
trum of a system of size L by a sum of peaks, enumerated 
by a counting index a, with position fi Q , width fj a , weight 
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Fig. 6. (Color online) Comparison of CheMPS, Bethe Ansatz 
and tDMRG + prediction, for HAFM structure factors with 
(a) k = 7r and (b) k = n/2. Dashed lines: Bethe Ansatz re- 
sults for L — oo, from Ref. HE Dashed-dotted lines: tDMRG 
results, from Ref. I22I Other lines show CheMPS results for 
L — 50, 100, 200, 300, reconstructed using a fixed ratio of 
N/L, namely 0.42 for (a) and 0.67 for (b). Circles mark 
Chebyshev peak maxima, also for L — 66, 150 and 250, for 
which no curves are shown. The lower inset in (a) zooms into 
the peak region using a linear scale, illustrating overbroad- 
ening. The upper insets of (a,b) show the peak heights vs. 
L (circles), together with a fit to the Bethe- Ansatz expecta- 
tion from Eq. (|46[) (dashed line) or to a straight line (solid 
line). In (b), triangles mark the lower threshold frequencies 
for which S(tv/2, uj) equals a fixed, small value, arbitrarily cho- 
sen as 0.1. Their 1/L — > extrapolation, shown in lower inset, 
gives an estimate for the lower threshold frequency, namely 
uii/n = 0.496 (marked by a square); the exact value is 1/2. 

W a and Gaussian shape K J (cf. Eq. (|18a[) ) : 

S £ (fc,w)^£W a i?£(w-fi a ) ■ (49) 

a 

Its Chebyshev reconstruction with Jackson damping, say 
Sft(k, oj), will have the same form, except that the peaks 
will be broadened to have widths, r\ a = (?y 2 + ^q) 1 ^ 2 , 



as explained before Eq. ([T9]) . If N is large enough, the 
broadened peaks will still be clearly separated (as for 
N = 999 or 250 in Fig. By fitting each peak (sep- 
arately, one by one) to a Gaussian, one can determine 
its position Q a , weight W a and effective width rj a , and 
deduce the intrinsic width via fj a = (rfc — rfa ,-J 1 / 2 . Wc 
find (not shown) that the intrinsic width grows with in- 
creasing frequency fl a . This implies, not unexpectedly, 
that higher-lying spin-wave excitations have shorter life- 
times. However, it also implies that higher-lying peaks 
eventually start to overlap, so that the analysis to be 
described below is feasible only for a limited number of 
low-lying peaks. 

The discrete peaks suggest a natural partitioning of 
the frequency spectrum into intervals I a : each con- 
tains one peak of weight W a at position £l a , extends 
halfway to the next peaks at fl a ±i on either side, and 
has width A Q = (fl a +i — f2 Q _i)/2. The first inter- 
val above the lower spectral threshold (wi) is defined 
slightly differently: 1\ has lower bound lo\ and width 
Ai = (a 1 + fl 2 )/2-u} 1 . 

Now, to produce a smooth curve devoid of finite-size ef- 
fects, the subpeaks must be broadened until they overlap 
substantially. However, if the weights in two neighbor- 
ing intervals differ, say W a > W a +i, such broadening 
inevitably transfers weight from interval I a to I a +i, re- 
sulting in overbroadening. 

Such overbroadening can be avoided by construct- 
ing a discrete representation of the spectral function, 
<Sdis(fc, ^q), defined by the set of coordinates 

{(fi aj S a )}, with S a = S dis (k, n a ) = W a /A a . (50) 

The identification of S a with W a /A a follows from apply- 
ing the definition of a spectral function, namely spectral 
weight per unit frequency interval, to the interval I a . 

Figure [7] shows the resulting discrete data points for 
four different system sizes. Remarkably, they all fall onto 
the same curve, which agrees well with the Bethe Ansatz 
result (dashed line). In particular, the first two or three 
data points for each L lie right on top of the Bethe Ansatz 
curve (dashed), see Fig. [7l left inset, beautifully mapping 
out the true shape of the spectral function down to the 
lowest discrete excitation frequency Cl a that exists for 
that L. Evidently, the discrete spectral function is com- 
pletely free from broadening artifacts, in marked contrast 
to the optimally broadened curves shown for L = 200 and 
400 (solid lines) (compare also Fig.^f)). This advantage 
comes at the price of specifying the spectral function only 
at discrete points, not via a continuous curve. However, 
for a system of finite size, such discreteness is fundamen- 
tally unavoidable. The good news is that the continuum 
curve S^(k, w) is evidently well mimicked by the discrete 
representation {(Oa, S a )}, and that CheMPS allows the 
latter to be determined in a straightforward fashion for 
system sizes well beyond what can be done with exact di- 
agonalization. We are not aware of any other numerical 
many-body method capable of doing so for system sizes 
as large as those considered here. 
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Fig. 7. (Color online) Discrete representation (Eq. (|50[0 of 
the structure factor S a = <Sdis(vr, Q a ) (symbols), for five dif- 
ferent system sizes. (The lower panel uses an enlarged verti- 
cal scale, to zoom in on the tail region.) For comparison, the 
Bethe Ansatz result (dashed line) and two optimally broad- 
ened spectra, for L — 400 and 200 (solid lines), are also shown. 
Left inset: zoom to low frequencies, showing that the dis- 
crete data completely avoids overbroadening, in contrast to 
the optimally broadened spectra. Right inset: log-log ver- 
sion of main plot. The frequency range does not extend low 
enough to be able to uncover the pure asymptotic predicted 
by Eq. ([46b]) , 



For larger frequencies the scatter of the discrete data 
w.r.t. the Bethe Ansatz curve increases, reflecting the fact 
that subpeaks begin to overlap there, making the extrac- 
tion of discrete data increasingly difficult. However, this 
is not a serious concern, since in this frequency regime 
optimal broadening is able to produce smooth spectra in 
good agreement with Bethe Ansatz anyway. 

To conclude this subsection, let us summarize the two 
main results of our finite-size analysis. The first con- 
cerns physics: for a chain of finite chain of L sites, the 
structure factor is dominated by a set of O(L) sharp sub- 
peaks, whose spacing and weight scale as 1/L. The sec- 
ond concerns methodology: CheMPS very conveniently 
allows this structure to be revealed or hidden, by sim- 
ply choosing N appropriately. Moreover, it can exploit 
information on the positions and weights of the discrete 
subpeaks to largely eliminate broadening artefacts. 



E. Comparison of tCheMPS to tDMRG 

Another possible scheme for reducing finite-size effects 
is to work in the time domain using linear prediction, as 
shown in Ref. for the HAFM. The idea is to calculate 



the Fourier transform of S(k,u>), namely 

L 

S(k, t) = J2 e lkix ~ X ' ] (S* (0)> , (51) 

x=l 

with x' chosen near the middle of the chain, and t chosen 
small enough that the spin excitation created at x' does 
not reach the edge of the system within t. The function 
S(k, t) thus obtained will contain only weak finite-size 
effects. It is then extrapolated to larger times via linear 
prediction techniques^ - — exploiting the fact that mo- 
mentum excitations typically exhibit damped harmonic 
dynamics, whose time-dependence can be extrapolated 
quite accurately. Since the extrapolated function extends 
to very large times, its Fourier transform yields good 
spectral resolution at low frequencies^ (with an accuracy 
that depends on that achieved during linear prediction). 

In Ref. the input correlator needed for linear pre- 
diction, S(k,t), was calculated using tDMRG. (Two 
examples of the resulting spectra arc included in our 
Fig. [6l) We note that S(k,t) can also be calculated us- 
ing CheMPS in the time-domain, to be called tCheMPS. 
Indeed the numerical cost for calculating S(k, t) by evalu- 
ating the requisite correlators (S x (t)S x >(0)} via Eq. ([28]) 
is essentially the same as calculating its Fourier transform 
S(k,u>) via Eq. ([23]) . since the corresponding Chebyshev 
moments \i n can be calculated using the same recursion 
scheme. In fact, if one defines Sk in Eq. (|43|) using a pure 
exponential e lh i instead of a sin function, the Chebyshev 
moments needed for S(k, t) are simply linear combina- 
tions of those of S(k, a;). 

To gauge the accuracy of tCheMPS, we have calculated 
S(n/2,t) using both tCheMPS and tDMRG. Figure^a) 
compares the results, and Fig. [8]Jb) characterizes the dif- 
ferences. We view the tDMRG results as benchmark, 
because for the times of interest, we have checked them 
to be well converged (with errors < 10~ 3 for t < 50, see 
Fig- IHb), dashed-dotted line). As expected, the agree- 
ment between tCheMPS and tDMRG is better for larger 
m. The differences are very small, but grow with time, 
from being (for m = 64) below 10 -3 for t < 10 to around 
10~ 2 for t ~ 30, beyond which finite-size effects start to 
appear. 

More generally, the results of Fig. [5] illustrate that 
CheMPS offers a viable route to time evolution for situ- 
ations where extreme accuracy is not required. Further 
comments on this prospect are included in the outlook, 
Sec. EE 



V. ERROR ANALYSIS 

The convergence properties of a Chebyshev expansion 
are mathematically well controlled and understood (see 
Eq. (fl4)) ). provided that the Chebyshev moments fi n are 
known precisely. Their evaluation via CheMPS, however, 
introduces various sources of numerical errors. This sec- 
tion is devoted to an analysis of these errors. In par- 
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Fig. 8. (Color online) (a) Time dependence of S(tt/2, i), calcu- 
lated with tCheMPS (lines) and tDMRG (symbols). Solid and 
dashed lines show, respectively, the real and imaginary parts 
of S. (b) The differences between tCheMPS and tDMRG 
(with a specified truncation error of 1(T 6 ) of S(k,t) for two 
values of k, and two values of m (dashed/solid). To estimate 
the accuracy of tDMRG, we also show (dashed-dotted) the 
differences between two tDMRG calculations performed with 
different truncation error thresholds, namely 10~ 5 and 10 -6 , 
requiring up to m = 75 or 125 states, respectively. 



ticular, wc seek to determine appropriate choices for the 
control parameters associated with the various CheMPS 
tasks listed in Table [U We perform this analysis mostly 
for a resonant level model (RLM), describing three local 
levels coupled to a fcrmionic bath. This model is intro- 
duced and discussed in App. El which, for the sake of 
completeness, also includes CheMPS expansions of the 
corresponding spectral functions. However, the details 
presented there are not needed for the following discus- 
sion. 

For the RLM, on the one hand, the CheMPS evalua- 
tion of the /z n is feasible to arbitrarily high orders, and on 
the other, exact diagonalization (to be denoted by sub- or 
superscript ED) of the single-particle Hamiltonian allows 
both the spectral function and the Chebyshev moments 
|U„ to be found exactly. We use the RLM-parameters 
specified in App. [A] throughout and focus mainly on the 
properties of one of its correlators, A^-y (without display- 
ing corresponding sub- and superscripts), which is de- 
fined in Eq. ()A2|) and whose behavior is representative 
for that of A„. 



A. Definition of error measures 

We will analyse both [i n - and /i n -moments, calculated 
from Eqs. (|25j) and (j34|) . respectively. The differences 
between CheMPS and ED can be quantified by the error 
measures 



5 ED 



I CheMPS 
I -CheMPS 



« ED I 
u ED l 



n < N max , 
n < 27V max . 



(52a) 
(52b) 



Moreover, to characterize the accuracy of CheMPS mo- 
ments without referring to exact results, we also consider 



C 1 - I An - f^n | , n < N max . (52c) 
We will also use cumulative versions of these, namely 



, ED 
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(53a) 
(53b) 
(53c) 



Furthermore, we also introduce an integrated error mea- 
sure for undamped spectral functions (using Jackson 
damping would yield qualitatively similar error mea- 
sures): 



J du \A 2N ^{±uj) - A 00 {±oj)Y 



(54) 



Here we use ± for A ± {lS) spectra proportional to 8(±oj) 
(see Eq. (|A2|> ) ■ and employ /%-moments for n < N max 
and /ijj-moments for N max < n < 2N max during spec- 
tral reconstruction. (Note that A ED of Eq. (|53T>[) was 
constructed to reflect this combination of /i n and fl n .) 



B. Comparison of CheMPS and ED moments 

Figure [HI contains the results of our comparison of 
CheMPS and ED moments for a fixed set of CheMPS 
parameters, stated in the figure legend. Figures [9ja,b) 
show Chebyshev moments [i n and jl n , Figs. G2c,d) the 
n-dependent error measures, <5 ED , i5 ED and From 
Fig. [9jc) we note several points: (i) For n < iV max , 
the //^-moments from CheMPS and ED agree to within 
about 1%; this illustrates that CheMPS is able to gen- 
erate rather accurate results for several hundered mo- 
ments at modest computational costs, (ii) /^-moments 
are more accurate than /i„-moments; the reason is that 
each ^„-moment depends on only one Chebyshev vec- 
tor, whereas each /Et n -moment depends on two. (Note, 
though, that if spectral reconstruction is performed by 
employing both ^„-moments for n < N max and jl n - 
moments for n > N max (as done, e.g., for Figs. [5] and 
[T3| , the reduced accuracy of the /2„-moments is offset to 
some extent if damping factors g n are employed, since 
these decay to as n approaches N, see inset of Fig. [3]) 
(iii) The error measures <5^ h and (5 ED are of comparable 
magnitude; this implies that 5^ h is a useful error quan- 
tifyer if exact results are not available. 
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Fig. 9. Comparison of CheMPS and ED results for Chebyshev moments of the RLM spectral function A lx . (a,b) show ji n - 
and fl„ moments (Eqs. (I25|l . (|34|l ) and (c,d) the n-dependent error measures <5 ED , 5™ and (Eqs. [52]) . plotted in (a,c) for 
n < -/V max = 200 and in (b,d) for JV max < n < 2iV max . In (b), the increase in moment magnitude starting around n ~ 250 marks 
the onset of resolving finite-size structure in the spectral function. (e,f) show the cumulative error measures A ED , A ED , A Ch 
(Eqs. (53])) and A A (Eq. (54])) for various combinations of the MPS dimension m, the number of energy truncation sweeps ns 
and the Krylov subspace dimension dx. 



The way in which theses errors depend on the various 
CheMPS control parameters can conveniently be anal- 
ysed using the cumulative error measures A ED , A ED , 
A Ch and A A . These are shown in Figs. [9^e,f ) for various 
combinations of m, ns and dx ■ Several observations can 
be made: (iv) When increasing the Krylov subspace di- 
mension dx , all cumulative errors decrease from dx = 20 
to 30, but the decrease saturates beyond dx = 30. (v) 
Increasing the number of energy truncation sweeps be- 
yond ns = 10 does not necessarily reduce the cumulative 
errors; on the contrary, most actually increase, implying 
that energy truncation sweeping should not be overdone, 
(iv) The cumulative errors depend only weakly on the 
MPS dimension m (except for dx = 10, which is unreli- 
able anyway), and tend to be smaller (!) for m = 32 than 
64 (compare Figs.^e) andOJf)). This trend suggest that 
the errors introduced by energy truncation grow if the 
mismatch between m and dx grows. Points (iv) to (vi) 
indicate that energy truncation is the limiting factor for 
reducing CheMPS errors, a fact that will be elaborated 
on in Sec. IV Cl below. 

To identify an optimal combination of CheMPS con- 
trol parameters, we have collected error data such as 
those shown in Figs. [9je,f) for each possible combina- 
tion of W* = (1.1, 1.5, 2.0), e t = (0.1, 0.01, 0.025), d K = 
(10, 20, 30, 50), ns = (5, 10, 20), and several m-values, for 
fixed maximum recursion number -/V max = 50 and con- 



vergence threshold A c = 10 8 . We concluded that the 
choices d K = 30, n s = 10, W* = 2W A and e t = 0.025 
robustly yield good results (also for the HAFM), and 
hence list these as recommended values in Table HI Ac- 
tually, the precise choice of et has only small effects on 
the error, as long as W, is chosen big enough. If W, is 
too small, however, the resulting spectral function will 
loose some weight at high frequencies, because numeri- 
cal errors may cause energy truncation to effectively also 
project out some contributions with energies smaller than 
the energy truncation threshold ep. 



C. Errors induced by recursion fitting and energy 
truncation 

To better understand the error dependence on m. dx 
and ns observed in points (iv) to (vi) of Sec. IV Bl above, 
let us analyse in more detail the errors generaturcd dur- 
ing recurrence fitting (Sec. ImA"]) and energy truncation 
(Sec. IIII Bp . The error incurred when constructing \t n ) 
from and |t n _2) using recurrence fitting is char- 

acterized by the relative fitting error Aj~ t = Afi t /|| \t n }\\ 
(Eq. ([31"]) ). The effect of projecting out high-energy states 
using energy truncation, \t n ) h-> Ptrl^n), can be charac- 
terized by the average truncated weight per site during 
one truncation sweep, N™ eep (Eq. (|4H]) ). and by the rel- 



16 



ative truncation-induced state change A£ r = A tr /|||i„)|| 
(Eq. (|4Tj) ). The latter measures intended changes in the 
state due to the truncation of high energy weight, but 
also incorporates the effects of unavoidable numerical er- 
rors. 

These quantities are analysed in Fig. [ljj]in dependence 
on to, d,K and ns- Continuing our list of observations 
from the previous subsection, we note the following fea- 
tures: (vii) Both Ag t and AJ; r are smaller than 1% al- 
ready for to = 32 (Fig. 110( a)). in accord with similar er- 
ror margins for <5™ in Fig.(9^c). (viii) Both A| t and Aj. r 
decrease with increasing to, but A\ r does so more slowly, 
and its decrease seems to saturate beyond to = 64. This 
implies that energy truncation is the main limiting fac- 
tor for CheMPS. The reason is that the intended purpose 
of energy truncation, namely to strip \t n ) from its high- 
energy components, modifies it in a way whose errors 
cannot be reduced to arbitrarily small values. Indeed, 
this is illustrated by the following two points: (ix) While 
both Atr and N™ eep initially decrease with increasing 
Krylov subspace dimension da, the decrease saturates 
for d K > 30 (Fig. Mb,c)); (x) While iV* r WGOp initially 
decreases with the number of sweeps ns, the decrease 
saturates already for ns i5 10 Fig. [TOTc). Qualitatively, 
the behavior shown in Fig. flUTc) is robust. (However, 
the choices of other CheMPS control parameters do in- 
fluence its quantitative details, such as the dx beyond 
which N™ ccp becomes (^-independent.) The lack of 
saturation of N^ ccp with ns implies that there is no 
automatic stopping criterion for truncation sweeps. In- 
stead, the choice of ns can be optimized as described 
in Sec. IV Bl where we already concluded that taking ns 
much larger than 10 actually deteriorates the results. 

Of course, truncation-induced errors can be avoided by 
simply using the full bandwidth, M4 = W, for which no 
trunctation is necessary. However, in our experience the 
gain in resolution obtained by using, instead, an effective 
bandwidth W* <C W, outweighs the small loss in accu- 
racy incurred by the necessity to then perform energy 
truncation. 

VI. DENSITY MATRIX SPECTRA 

The effects of energy truncation can be understood in 
more detail by considering the reduced density matrix 

p n = Tr ha lf I ( 55 ) 

where the trace is over one half of the chain. Let us anal- 
yse the n-dependence of the spectrum of its eigenvalues, 
say p n (i). It can be used to quantify the entanglement 
encoded in \t n ), via the associated entanglement or bond 
entropy, 

C nd = -E^)H^))- (56) 

i 

Figure [TT1 shows such density matrix spectra for both the 
RLM (panels (a,b)) and the HAFM (panels (c,d)), calcu- 
lated using both the full many-body bandwidth W* = W 



(panels (a,c)) and a smaller effective bandwidth W, (pan- 
els (b,d)). The n = line in all panels shows the eigen- 
value spectrum po{i), which reflects the entanglement en- 
coded in \t ) = C\0) at the start of the recursion pro- 
cedure. In principle one would expect the entire spec- 
trum of density matrix eigenvalues p n (i) to shift or rise 
to higher values as n increases, since multiplying |t n _i) 
by H' when calculating \t n ) (cf. Eq. (|30|) ) generates en- 
tanglement entropy. Such a spectral rise with increasing 
n is indeed observed in all four panels of [Til but the rise 
eventually saturates for sufficiently large n. The speed of 
the initial stages of the rise differs from panel to panel. 
For the density matrix spectra calculated without energy 
truncation (Fig. Illf a.c)), the initial rise is rather slow, 
in particular for the RLM (Fig. [TTl a). where the rise 
is preceded by a slight initial decrease), reflecting the 
lack of strong correlations of this model. In contrast, for 
density matrix spectra calculated with energy truncation 
(Fig. Illf b.d)). the initial rise is very rapid, and its sub- 
sequent saturation sets in at quite small n (of order 20 
to 30). Thus, energy truncation evidently has the effect 
of increasing entanglement entropy. The reason is that 
the latter is calculated in a different basis (the eigenbasis 
of p n ) than that used to perform energy truncation (the 
local eigenbasis of H'). 

According to Fig. [TJJd), the small MPS dimension of 
m = 32 used for the HAFM in Fig.[S{a) in effect amounts 
to discarding the contributions to the reduced density 
matrix of all states with weight below a threshold of 
around 10~ 3 . This threshold is rather large compared to 
typical DMRG calculations, where characteristic trunca- 
tion errors lie in the range 10 -6 to 10 -8 . It is remarkable 
that CheMPS is nevertheless able to give rather accurate 
results (such as reproducing CV results obtained using 
to cv = 1500). 

This efficiency appears to be an intrinsic feature of 
CheMPS, arising from the recursive manner in which the 
Chcbyshev vectors \t n ) are constructed. Evidence for this 
conclusion is presented in Fig. [T2T a'). which shows the 
bond entropy S'^' ond associated with \t n ) as function of re- 
cursion number n. Remarkably, the bond entropy shows 
no tendencies towards unbounded growth, even up to val- 
ues as large asn = 2000. Quite to the contrary: although 
the bond entropy increases somewhat when increasing m 
from 32 to 128 (with W, = 6.3), for either case it tends 
to decrease with recursion number n, and similarly for 
the choice W* = W without energy truncation. All of 
this is very encouraging, since it indicates that n can be 
increased, apparently at will, without incurring any run- 
away growth of DMRG truncation errors. The reasons 
for this fact will be recapitulated in the summary below. 

For comparison, Fig. I12f b) shows the bond entropy 
St ond of a tDMRG calculation of the time evolution of 

\tp(t)) = e~ lHt S x= 5o\0}- This entropy is, overall, smaller 
than the 5'^ ond of the Chebyshev vectors, because the ini- 
tial state for the time evolution involves an excitation at 
only one site, whereas the starting state for the CheMPS 
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Fig. 10. (Color online) (a) Relative fitting error AJi t = Aat/ |||tn)|| (Eq. (J3TJ) an d relative truncation-induced state change 
A[ r = Atr/||jin)|| (Eq. (J4TJ), as functions of recursion number n, for three different choices of MPS dimension m. Both 
quantities decrease with increasing m, but more strongly so since recurrence fitting is a strictly variational procedure, 
whereas energy truncation is not. (b) A£ r as function of n, and (c) the average truncated weight per site N^ eep (for n = 20) 
as function of truncation sweep number k. Both (b) and (c) show results for four choices of Krylov subspace dimension (Ik, 
whose dx-dependence saturates beyond cLk = 30. 
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Fig. 11. (Color online) Eigenvalue spectra p n (i) of the re- 
duced density matrix at the center of the system for several 
expansion vectors \t n ) of (a,b) the RLM with Lt = 101, and 
(c,d) the HAFM with L = 100. In (a,c) we used the full many- 
body bandwidth W, — W without energy truncation, in (b,d) 
a reduced effective bandwidth with energy truncation. 



recursion involved a linear combination of local excita- 
tions, Sk\0) (see Eq. ([43f ). The most striking difference 
between S^ ond and S^° nd , however, is that the former 
shows no trend to increase with n, whereas the latter 
does with t. The increase in S^° nd occurs in spurts, that 
happen each time a spin wave gets reflected from one of 
the ends of the system, at which point more numerical re- 
sources are required to keep track of the superposition of 
incident and reflected spin waves. For the present prob- 
lem, the increase in S^° nd was not severe and remained 
completely under control (staying below S^ oad through- 
out). Nevertheless, we do believe that the contrast be- 
tween Fig. [T27 a) and Fig. H^T b). showing a nonincreasing 



trend for 5^ ond vs. an increasing trend for S^ ond , is strik- 
ing and significant. It suggests that for situations that 
feature strong entanglement growth with time, tChcMPS 
might be a promising alternative to tDMRG. 



VII. SUMMARY 

In this work, we have described CheMPS as a method 
for calculating zero-temperature spectral functions of 
one-dimensional quantum lattice models using a combi- 
nation of a Chebyshev expansion and MPS technology. 
To summarize our analysis, we would like to highlight 
what we believe to be the two most important features 
of CheMPS, namely its efficiency and its control of spec- 
tral resolution. 

Efficiency. The first main feature is that CheMPS pro- 
vides an attractive compromise between accuracy and ef- 
ficiency. It is capable of reproducing correction vector 
results in the frequency domain and tDMRG results in 
the time domain with comparably modest numerical re- 
sources. In particular, surpringly small values for the 
MPS dimension of m are sufficient, even for obtaining 
spectral resolution high enough to resolve finite size ef- 
fects in great detail. (For example, m = 32 sufficed for 
the spin-i antiferromagnetic Hcisenbcrg model.) This re- 
markable efficiency, which we had not anticipated when 
commencing this study, appears to be a consequence of 
several factors: (i) CheMPS does not suffer from a run- 
away growth of DMRG truncation error with increasing 
n, because the information needed to construct the spec- 
tral function with a specified accuracy, say 0(1/N), is 
not encoded in a single state, but uniformly distributed 
over N distinct Chebyshev vectors \t n ). (ii) These can 
be determined from Chebychev recurrence relations in- 
volving only three terms, so that it is never necessary 
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Fig. 12. (Color online) Entanglement or bond entropy S hond 
for the k = tt spectral function of the HAFM. (a) S^ ond for 
the Chebyshev vectors \t n ) and (b) S t bond during the tDMRG 

time evolution of e _lHt S , a;= 5o|0). In both (a) and (b), solid 
and dashed lines show the maximum bond entropy and the 
bond entropy at the middle of the system, respectively, (a) 
iSn is shown for two choices of 114; the dotted line is from 
a calculation with a reduced m — 32 and some entropy is 
lost due to truncation. The red arrow marks the expansion 
order roughly necessary to reach the time t — 85 using the 
tCheMPS technique for W* = 6.3, here N t < SB = 271. To reach 
the same time using W* = W = 68.8 an order of expansion of 
Nt<gs = 2961 would be necessary, (b) S^° nd is shown for two 
choices of the truncation error e p . 



to accurately represent the sum of more than two MPS. 
(iii) Moreover, these recurrence relations are numerically 
stable, i.e. the inaccuracies in the calculation of Cheby- 
shev vectors \t n ) do not cause the Chebyshev expansion 
to diverge, (iv) Finally, the accuracy needed for each \t n ) 
is set by that needed for fi n = (0|B|t„) (Eq. ([25])). which 
does not need to be better than the specified accuracy, 
namely (D(l/N). 

For spectral functions with a finite spectral width 
Wa (which is typically much smaller than the many- 
body bandwidth W), ChcMPS offers a further attrac- 
tive feature for enhancing efficiency: one may use an "ef- 
fective bandwidth" W* of order Wa (we typically take 



W, = 2W/\), which enhances spectral resolution by a 
factor W/W,, at the cost of requiring additional energy 
trunctation sweeps. The latter are not necessary if one 
takes W* = W, but then considerably higher expansion 
orders are necessary to achieve comparable resolution. In 
our experience the benefits of enhanced resolution offered 
by the choice W* = 2Wa outweigh the costs of energy 
truncation. 

Control of spectral resolution. The second main feature 
of CheMPS is that it offers very convenient control of the 
accuracy and resolution of the resulting spectral func- 
tion, by simply adjusting the expansion order TV. This is 
particularly useful for studying finite-size effects, as ex- 
emplified in Fig. 03 On the one hand, Fig. E^b) shows 
very strikingly that the structure factor an HAFM chain 
of finite length is dominated by a set of discrete subpeaks 
which may be associated with the quantized eigenener- 
gies of spin wave excitations in a finite system. CheMPS 
allows the energies and weights of these excitations, and 
their dependence on L, to be determined with unprece- 
dented accuracy and ease, by simply increasing TV until 
the peaks are well resolved. On the other hand, Fig. [5jf ) 
shows that the limit L — > oo may be mimicked by choos- 
ing N just small enough that the finite-size subpeaks are 
smeared out. Though the peak shape thus obtained is 
slightly overbroadened (see inset of Fig. [5ff )), this over- 
broadening can be eliminated completely (see Fig. [7]) by 
using a discrete representation of the spectral function, 
that uses the energies and weights of the discrete sub- 
peaks as input. The ability to fully eliminate overbroad- 
cning effects even for very large many-body systems is, to 
the best of our knowledge, a unique feature of CheMPS. 

On a technical level, the implementation of CheMPS 
requires only standard MPS techniques, such as the ad- 
dition of different states and the multiplication of opera- 
tors. For energy truncation, single-site sweeping needs to 
be set up with a new kind of local update, as described 
in Sec. lIIIBl However, this procedure is not too different 
from other known local update prescriptions and can be 
implemented with modest programming effort. 



VIII. OUTLOOK 

Regarding future applications of CheMPS, two direc- 
tions for further methodological development appear par- 
ticularly promising, namely time-dependence and finite 
temperature. A few comments are due about each. 

Time dependence. While the good agreement between 
tCheMPS and tDMRG reported in Fig.[5]is encouraging, 
a detailed analysis of tCheMPS should be performed to 
understand the nature of its error growth with time, and 
to explore under which conditions, if any, tCheMPS of- 
fers competitive advantages relative to tDMRG. On the 
one hand, tDMRG has the advantage that highly effi- 
cient Krylov methods can be used to optimize the evalu- 
ation of e~ %H \ip{t)) w.r.t. the state \4>{t)) being prop- 
agated; however its numerical costs increase rapidly if 
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\tj}{t)) contains a broad spectrum of excited states. On 
the other hand, CheMPS has the advantage (i) that the 
Chebychev expansion of the operator e~ lHt can be ap- 
plied with equal accuracy to every state in the Hilbert 
space, in particular also highly excited ones. Moreover, 
(ii) very large evolution times might be achieved more 
easily with tChcMPS than tDMRG, since the former rep- 
resents \ip(t)) as a sum over many Chebyshev vectors (see 
Eq. dU)), thereby being potentially less susceptible than 
tDMRG to the growth of truncation errors (as discussed 
in the introduction, and exemplified in Fig, \12\i . We ex- 
pect that for some applications (i) and/or (ii) may offer 
advantages for tChcMPS over tDMRG, e.g. for calculat- 
ing quantum quenches starting from strongly noncquilib- 
rium initial states, but leave a detailed investigation to 
the future. 

Finite temperature. The fact that CheMPS uniformly 
resolves the entire energy spectrum of H suggests that it 
should be particularly suited for calculating the spectral 
functions A% c (uj) = J ^e lut G% c (t) of finite temperature 
correlators such as 

G BC (t) = Tr[p T B(t)C(0)], Pt = —- (57) 

According to Ref. [f| such a spectral function can be eval- 
uated using Chebyshev expansions by proceeding as fol- 
lows: Express the partition function as 

Z = J due-P" p(u)) , (58a) 

by introducting the density of states 

p(u) = Tr[<y(w - H)} , (58b) 
and the spectral function as 

A% c {uj) = ^ J dCoe- 0Q p BC {uj,uj + uj) , (59a) 

by introducing the density of matrix element s 39 ^ 

p BC (Cj,Lj) = Tr[5(Cj-H)B5(uj-H)C} . (59b) 

Then Chebyshev expand the 6- functions in Eqs. (|58b|) 
and (|59b|) using Eq. ([25)) (after suitably rescaling Hamil- 
tonian and frequencies) . The resulting Chebyshev expan- 
sions will contain moments of the form 

fx? = Ti[T n (H')} , (60a) 
fJ-nn' = Tr[T n (H')BT n ,(H')C) . (60b) 

We now note that this framework is very well suited for 
an MPO implementation, which would consist of three 
steps: (i) Using Chebyshev recurrence relations, recur- 
sively construct and store MPO representations for each 
operator T n (H'); we expect (based on our experience 
with the Chebyshev vectors \t n )) that this should be pos- 
sible without runaway costs in numerical resources, since 
the construction of T n (H') requires only H'T n _i(H') and 
T n _2(H'). (ii) Calculate the moments in Eqs. (|6T)|) by 



evaluating the traces, which is straightforward in the con- 
text of MPS/MPO. (iii) Insert the resulting moments into 
the reconstructed Chebychev expansions for p(uj) and 
p BC (ui,ui), and finally evaluate the integrals Eqs. (|58aD ) 
and (|59a|) . Note the economy of this scheme: after once 
constructing the MPO for each T n (H'), and once evalu- 
ating the trace for each moment pP n and , the spectral 
function A BC (u>) can be calculated for arbitrary combi- 
nations of oj and T. The implementation of this strategy 
is left for future studies. 

We conclude by remarking that the idea of using 
Chebyshev expansions in the context of many-body nu- 
merics, advocated in inspiring fashion in Ref. [H, can be 
implemented in combination with any method that is 
able to efficiently apply a Hamiltonian H to a state 
Chebyshev expansions optimize the resolution that can 
be extracted from a limited number of applications of 
H. While CheMPS is based on doing this using MPS 
methods for one-dimensional lattice models, similar de- 
velopments have been pursued within the context of ex- 
act diagonalizatio n 41 ' 42 and Monte Carlo^ methods, and 
Chebyshev expansions should also be useful in combina- 
tion with tensor network methods for two-dimensional 
quantum lattice models. 
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Appendix A: Resonant level model 

This appendix introduces the fcrmionic resonant level 
model (RLM) that was used for the error analysis of 
Scc.[V] and presents CheMPS results for its spectral func- 
tions. 

The RLM is defined by the following Hamiltonian: 

"d »d Lb 

h klm = e + E v ' E + h - c -j ( A1 ) 

(=1 4=1 fc=l 

L b 

+ E e *4 5 /e ■ 
k=l 



20 



It describes a set of rid discrete, "local" , non-interacting 
fcrmion levels with energies Ei, that hybridize with 
strengths Vi with a band of Li, 1) fcrmion levels with 
energies assumed uniformly spaced within the inter- 
val [— W6,Wb]. We choose W& = 1 as unit of energy 
throughout this section. We will parametrize the hy- 
bridization strengths Vi in terms of the associated level 
widths r t = ir^V?. 

The spectral function Aij(ui) = Afj(u)) + A~j(uj) has 
two contributions, 



Ar.( w ) = ^K(-a;), (A2) 

describing particle and hole excitations, which at T = 
are proportional to step functions 9(±lo) that vanish for 
u < or > 0, respectively. Since the RLM Hamiltonian is 
quadratic, the problem can be solved by diagonalizing the 
single-particle problem. In the continuum limit L b — > oo, 
this yields the following exact expression for the spectral 
function^ for \oj\ < D b = 1: 



too 



lim 3 

r/->0+ 7T 



IT] 



a( W )]- 



(A3) 



a,» - -Vrj^ 





LJ — D b 






— in^j 




uj + D b 





where T and A are matrices of dimension x n^. The 
Chebyshev moments (j, n for the finite system of length L 
can also be found exactly, by evaluating the expectation 
values Eq. (|24|) using the (numerically-determined) exact 
single-particle eigenstates of H. 

The Hamiltonian (|Alj) corresponds to a "star geome- 
try", since each local level couples to every band level. 
For the purposes of using CheMPS, however, it needs to 
be transformed to a "chain geometry" of the form 



i=l i=l 
L b -l 



(A4) 



This can be achieved^ using Lanczos tridiagonalization 
of the band part of the Hamiltonian, thereby determining 
the hopping coefficients Xg. 

Starting from Eq. (|A4|) . we have used CheMPS to cal- 
culate the diagonal components Ajj of the RLM spectral 
function for a model with = 3 local levels. In con- 
trast to Section [TV CI our interest here is not in analysing 
finite-size effects, but in determining how the CheMPS 
parameters need to be adjusted to recover the exact con- 
tinuum function A°° of Eq. (|A3[) . Thus, we purpose- 
fully chose a set of model parameters leading to three 
well-separated peaks of slightly different widths, taking 
Sj e {-0.5,0.1,0.6} and Tj e {0.04,0.06,0.03}, and 
chose the number of band levels L b = 101 large enough 
that the finite-size spacing wi ~ l/£& = 0.01 is some- 
what smaller than the smallest peak width, T^. By choos- 
ing the expansion order for each curve such that the effec- 
tive broadening lies in the window between the finite-size 



spacing and the intrinsic peak width, ujl < t)n < Tj, it 
should be possible to reveal the shape of A°° quite accu- 
rately without yet resolving finite-size subpeaks (though 
traces of the latter might show up for ^33, for which 
this window is small). To this end, we used the follow- 
ing criterion for choosing A when reconstructing A^: 
the effective broadening t]n was taken as large as possi- 
ble without lowering the peak height significantly below 
that of A^ (this corresponds to choosing r)N ^ Tj). 

The results of these calculations are summarized in 
Fig. [13J all spectra shown there were obtained by per- 
forming separate expansions for the positive and negative 
branches, A^{ui) (with one exception, noted below). 

Figures [T37 al-gl) were calculated using an effective 
bandwidth of Vt4 = 2.0 (with e t — 0.025) for each branch, 
corresponding to roughly twice the spectral width of each 
branch, which is of order of the single-particle band- 
width, Wa ^ W b = 1. For this choice, an MPS di- 
mension of merely m = 32 was found to suffice for accu- 
rate recurrence fitting. Figure [T37al) illustrates a num- 
ber of points: (i) By choosing 77 according to the above 
criterion of recovering the correct peak height, excellent 
agreement with the continuum limit A 00 of Eq. (|A3[) is 
obtained over most of the frequency range, (ii) This is the 
case both with and without Jackson damping (thin black 
or blue lines, respectively), but with Jackson damping, 
higher expansion orders are needed to obtain the cor- 
rect peak heights, since Jackson damping induces some 
artificial broadening (by a factor of n, sec Eq. (|18ap ). 
(iii) Small oscillations remain in some frequency ranges 
(see Figs. fT3T bl-gl) for zooms). These stem from three 
sources: finite-size subpeaks, numerical inaccuracies and 
step function artefacts near uj = (cf. points (iv), (vi) 
and (viii) below, respectively), (iv) For the spectrum 
with the narrowest peak, ^33, the window between ujl 
and T33 is so small that the criterion of reproducing the 
continuum peak height implies that small finite-size sub- 
peak remain visible, see Figs. [T37cl-gl) for zooms, (v) 
In contrast, such oscillations are almost entirely absent 
for the broadest peak, A22 (see Figs.n3Tdl.el)). since its 
width T2 is somewhat larger than w^. 

In order to illustrate the effect of energy truncation, 
Figs. [13fa2-g2) show the same spectral functions as 
Fig. HSJal-gl), but now setting W* = W, the full many- 
body bandwidth (here = 52.3), so that no energy trun- 
cation is needed. This allows us to make some additional 
instructive observations: (vi) Using the full bandwidth 
yields results of higher quality, in that numerical artefacts 
are significantly weaker (except near uj = 0), compare 
Figs. [T37d2-g2) and (dl-gl). The reason is that energy 
truncation constitutes CheMPS' s dominant source of er- 
ror (as shown in Section [V] below); its avoidance thus 
yields more precise Chebyshev moments [i n , especially 
for n > A^max. (vii) However, this improvement is nu- 
merically expensive: the increased effective bandwidth 
neccessitates larger expansion orders A, which in turn 
requires a higher MPS dimension (here m = 128). (viii) 
For the present model, it was possible to calculate sev- 
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L b =1 01 , e,=0.025, A c =10"' 
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Fig. 13. (Color online) Diagonal spectral functions Ajj(u>) of a 3-level RLM. Thick solid lines show the continuum limit, 
Ajj(u), from Eq. (|A3|I . Dashed and dashed-double dotted lines show CheMPS results, Ajjj or Ajj, with or without Jackson 
broadening, respectively, calculated for Lb = 101 band levels. For each spectrum, the effective broadening ryjv was taken as 
large as possible without lowering the peak height significantly below that of A°° . In (al-gl) we used an effective bandwidth 
of W* — 2.0 and in (a2-g2) the full many-body bandwidth W* = 52.3. The latter requires significantly larger expansion orders, 
but exhibits less numerical inaccuracies, compare (bl-gl) and (b2-g2), which represent zooms of the rectangles indicated in 
(al) and (a2), respectively. (b,c): Gibbs oscillations arise if A^ioj) are expanded separately, so that CheMPS attempts to 
resolve their 9{iiuj) steps. Expanding instead their sum, A + (ui) + A~ (— uj), and performing a Jackson-damped reconstruction, 
we obtain the smooth dashed-dotted line in (bl-el) (calculated using ED moments). (d,e) A22 nicely reproduces A22, because 
the peak is somewhat broader than ljl- (e-g) ^33 shows small but distinct finite-size wiggles, because the main peak is so 
sharp and narrow that recovering its height fully, requires an j/m so small that it is comparable to ujl- 



eral thousand moments without encountering numerical 
instabilities; this illustrates the fact that the Chebychev 
recurrence relations are numerically stable. 

Finally, let us address (ix) the wiggly artefacts near 
10 = 0. They reflect the fact that CheMPS was separately 
applied to the positive and negative branches of the spec- 
tral function, A ± (u>), shown in zooms in Figs. fT3T b) and 
Fig. [T37c). respectively. These are proportional to step 
functions 8(±uj), and hence abruptly dip to zero for uj < 
or > 0, respectively. The wiggly artefacts correspond 
Gibbs oscillations decorating these sharps dips. This 
problem can be avoided by performing a single Cheby- 
shev expansion of the sum, A + (uj) + A~ (—uj), which is a 
smooth function and leads to the perfectly smooth long- 
dashed line in Figs. I13f b.c). This improvement comes 
at roughly twice the numerical cost, since it requires a 
doubling of the spectral range to uj e [— Wa,Wa)- this 
implies a slight but obvious modification of the transfor- 



mations from uj to uj' and from H to H' to account for 
the shifted range of uj; and a doubling of W* and hence 
of the expansion order TV required to achieve a specified 
resolution. 

The main conclusions from our CheMPS calculations 
for the RLM are as follows: The strategy of using twice 
the spectral width as effective bandwidth (M4 = 2Wa) 
and performing energy truncation (Fig. I13f a)) is a satis- 
factory compromise between efficiency (only a few hun- 
dred Chebyshev moments arc needed) and accuracy (for 
which energy truncation is the main limiting factor). If 
desired, better results can be obtained by using the full 
bandwidth (144 = W) and thus avoiding energy trun- 
cation, albeit at the cost of significantly increasing the 
required expansion order by the factor W/2Wa- Never- 
theless, the calculation of Chebyshev moments /i n with 
very large n is feasible due to the remarkably numerical 
stability of Chebychev recurrence relations. 
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